import sympy as sp

def validate_cr_residuals(phi_expr, psi_expr, variables, test_points, extra_subs=None):
    """
    数值验证柯西-黎曼方程残差（支持极坐标/直角坐标系）
    
    参数:
        phi_expr (sp.Expr): 速度势表达式（实部）
        psi_expr (sp.Expr): 流函数表达式（虚部）
        variables (list): 变量列表，如 [x, y] 或 [r, theta]
        test_points (list): 测试点坐标，格式为 [(var1_val, var2_val)]
        extra_subs (dict): 额外符号替换（如 {a: 1.0, U: 1.0}）
    
    返回:
        dict: 各测试点的最大残差
    """
    residuals = {}
    var1, var2 = variables

    # 计算偏导数
    if variables[0].name == 'r':
        # 极坐标模式
        dphi_dr = sp.diff(phi_expr, var1)
        dpsi_dtheta = sp.diff(psi_expr, var2)
        cr1 = dphi_dr - (1/var1) * dpsi_dtheta
        dpsi_dr = sp.diff(psi_expr, var1)
        dphi_dtheta = sp.diff(phi_expr, var2)
        cr2 = dpsi_dr + (1/var1) * dphi_dtheta
    else:
        # 直角坐标模式
        dphi_dx = sp.diff(phi_expr, var1)
        dpsi_dy = sp.diff(psi_expr, var2)
        cr1 = dphi_dx - dpsi_dy
        dphi_dy = sp.diff(phi_expr, var2)
        dpsi_dx = sp.diff(psi_expr, var1)
        cr2 = dphi_dy + dpsi_dx

    # 数值评估
    max_error = 0.0
    for point in test_points:
        substitution = {
            var1: point[0],
            var2: point[1],
            **(extra_subs if extra_subs else {})  # 合并额外符号替换
        }
        # 强制转换为 Python 浮点数
        error_cr1 = float(abs(cr1.subs(substitution).evalf()))
        error_cr2 = float(abs(cr2.subs(substitution).evalf()))
        current_max = max(error_cr1, error_cr2)
        residuals[str(point)] = current_max
        if current_max > max_error:
            max_error = current_max

    residuals['max_error'] = max_error
    return residuals